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Abstract. The entanglement entropy of noninteracting fermionic system confined to 
two-dimensional (2D) honeycomb lattice on a torus are calculated. We find that the 
entanglement entropy can characterize Lifshitz phase transitions without a local order 
parameter. In the noncritical phase and critical phase with nodal Fermi surface, the 
entanglement entropy satisfies the area law. The leading subarea term is a constant 
in gapped phase rather than a logarithmic additive term in gapless regime. The 
tuning of chemical potential allows for a nonzero Fermi surface, whose variation along 
particular direction determines a logarithmic violation of the area law. We perform the 
scaling of entanglement entropy numerically and find agreement between the analytic 
and numerical results. Furthermore, we clearly show that entanglement spectrum is 
equivalent to edge spectrum. 
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1. Introduction 

Entanglement entropy is a unifying theme in the interdisciplinary theories during the 
past decade. It is a helpful bridge between quantum gravity, quantum information 
and condensed matter physics. It is defined as the von Neumann entropy of one of the 
reduced states for a pure state. Given a normalized wave function |H/ab) and a partition 
of the system into subsystems A and B , the associated reduced density matrix of A part 
is obtained from pa = Trs|'h J 4 s)(d / As| by tracing out the degrees of freedom in B part. 
The bipartite entanglement entropy is given by, 

S vN = -Tip A log 2 p A . (1) 

Related generalized quantities are the Renyi entropies defined as 

Sn = log 2 Trp^, (2) 

1 — 71 

whose limit recovers the von Neumann entropy via relation 5 v n = lim, w i S n . Once the 
complete values of S n for all n are known, one can obtain all eigenvalues of pa, which are 
called the entanglement spectrum. The spooky nonlocal nature makes the entanglement 
entropy hold great advantage in the context of condensed matter physics. On one hand, 
it plays an essential role in determining the efficiency of numerical algorithms based on 
matrix product states, as well as an extension to tensor network states [TJ and Projected 
Entangled Pair States (PEPS) [2j. On the other hand, it is proven that the entanglement 
entropy is able to capture the quantum criticality in a many-body Hamiltonian M- 
The constituents become extraordinarily entangled at criticality. Quantum critical 
points are successfully pinpointed in both free systems 0 0 and interacting systems 
00 0 EDI HU, especially for those critical ground states without prior knowledge of 
order parameters, e.g., disordered systems[T2j 131 Q3J and quantum spin liquids [ 15] . 
In particular, the discovery of a topological phase transition, where distinct phases 
of matter have the same symmetry, was in urgent need of basis-independent way of 
identifying quantum critical points. In this respect, extensive work has utilized the phase 
sensitivity of the entanglement entropy to successfully characterize the topological order. 
The universal constant term in the bipartite entanglement entropy has been shown to 
be a probe of topological contribution mm- The success of entanglement entropy 
as a diagnostic tool is essentially rooted in its peculiar scaling properties in a quantum 
state. The growing awareness of the monotonic scaling of von Neumann entropy can 
be traced back to the study of black hole physics [THl US |20j, and then the profound 
relations to conformal field theory (CFT) in (l+l)-dimensions were established [2T[ .22]. 
Despite the formal simplicity, it is surprisingly difficult to calculate the entanglement 
entropy analytically, even for noninteracting models. The entanglement entropy in one¬ 
dimensional system has been extensively studied and the scalings have been explicitly 
established by CFT [3J, [23]. Intensive research in terms of numerical algorithms unveiled 
a few qualitative and quantitative aspects in higher dimensions pHl 25] . 

In this work, we study the half-space entanglement entropy of free fermion in 
the anisotropic honeycomb lattice with fine tuning chemical potential. In order to 
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discriminate the relation between various boundaries and entanglement entropy, we 
place the fermion system on a torus. The scaling behavior of the entanglement entropy 
in a two-dimensional hexagonal lattice demands a comprehensive investigation. The 
basic interest in honeycomb lattices comes from the graphene and associated Dirac 
fermion physics. More interestingly, a recent influx of ideas and toolboxes can create 
the so-called ” artificial honeycomb lattice” [26], such as photonic lattices m, coupled 
molecular systems [28] as well as optical honeycomb lattice [29]. The manipulations of 
Dirac points are possible through these artificial methods j30[ 31,132]. The asymmetric 
tunneling amplitudes cause a topological Lifshitz transition from a semimetal to an 
insulator, which is characterized by the topological change of the Fermi surface in the 
Brillouin zone [33]. It was found that the Lifshitz points can be also captured by the 
entanglement entropy [34]. Moreover, based on the bulk-edge correspondence [35] [36], 
the entanglement entropy and the Berry phase are closely related in a topologically 
nontrivial system through the zero-energy edge states [37]. A fascinating aspect of 
graphene-like structure is that there are two types of edges, i.e., zigzag and armchair 
edges. The role of edges in determining the entanglement entropy in honeycomb lattice 
and associated scaling behavior are not fully understood. In addition, the chemical 
potential can be well tuned in graphene by a gate voltage, so the dependence of the 
entanglement entropy on the chemical potential is a nontrivial issue to be discussed. 


2. Tight binding model on honeycomb lattice 


The anisotropic fermionic tight-binding model on a two-dimensional honeycomb lattice 
reads 

H = - ^2 ^2 **(4Vh7! + bt+r^r) - h y" a\a r - n ^ b%, (3) 

t£Aa i=a,b,c tGAa rEA b 

where a r (b r ) is the annihilation operator on the triangular sublattice A a o (respectively 
A b •) of the rth dimer, r t is the nearest neighbor bond vector, and /i is the chemical 
potential. For simplicity, we confine the hexagonal lattice to the brick lattice without 
changing lattice topology as shown in Fig[U(a). The nearest-neighbor separations are 
given by 


Ta= (0,1), 7* =(-1,0), T c = (1, 0), (4) 

where the bond length is set to unit. Taking Fourier transformation for all lattice sites, 
the Hamiltonian ([3]) reduces to 


H = 



iA(jfe) 
-f(k) - n 


C k-> 


( 5 ) 


and then the eigenspectrum is 

e(k) = ±\J £ 2 (fc) + A 2 (k) - n, (6) 

where £(&)= t a cosk x + cos k x + t c cosk y , A(k)—t a smk x - 4sin k x -t c smk y , c g=( c p 
ct _) and k = (k x , k y ) is the Brillouin zone shown in FigJT](b). The bulk band structure 
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in Eg. (15]) is invariant under inversion : k —y —k. In the parameter space | t a — tf\ < t c < 



Figure 1. (Color online) (a) Schematic illustration of brick-type graphene torus. A 
virtual bipartition with zigzag edges is made horizontally. The shaded region denotes 
one of the subsystems. The dashed (red) box indicates a unit cell. A bipartition with 
armchair edges can be implemented by truncating the lattice vertically, (b) BZ sq is 
the Brillouin zone of square lattice, and the double count for k and k + tt indicates the 
BZ (region 1-4) is the Brillouin zone of the hexagonal lattice. 


t a + tb (containing symmetric case t a = h = t c ), the valence and conduction bands 
intersect linearly at two /cD-p°ints (known as the Dirac points), given by 

m2 _ m2 _ m2 \ 

' b r L n b U \ \ 

(7) 


kn = 


2n — arccos 


arccos 


2 t a t b > 

{tl+tptl-hl-tp 2 ■ 

2 tatbtZ 


which will be suited at Appoints, i.e., (±| 7 r, 0), of the first Brillouin zone for symmetric 
hoppings t\ = t 2 = £ 3 . The spectra and the corresponding Dirac points for a semimetallic 
phase is displayed in Figj2](a). The low-energy electronic states at any Dirac point 
are described by the massless relativistic Dirac fermions, which lead to a number of 
unconventional and fascinating phenomena, such as anomalous quantum Hall effect 
[38] . The asymmetric hoppings make the two Dirac points move away from the K- 
points and approach each other [3Hj. At a critical asymmetry the Dirac points merge at 
one of the four inequivalent M-points, that is, (± 7 r, 0) and ( 0 ,± 7 r). A band gap opens 
upon increasing the asymmetry further; see illustration in Figj2](b). The topological 
invariant that characterizes the topological properties is given by the Zak phase, which 
is the integration of Berry connection over in the reduced Brillouin zone [ 401 14T ] . If the 
integral loop encloses a gap closing point, the Zak phase gives 

If — — i f (v-(k)\d^v_(k)) ■ dk = tt, ( 8 ) 

J —TT 


where v-(k) = 


is the eigenvector of lower band and 9k = Arg[£(fc) +iA(k)}. 


Otherwise it is zero. The merging is a topological transition since the Berry phases ± 7 r 
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Figure 2. (Color online) Three-dimensional view of the energy spectrum of the 
Hamiltonian ([3]) for (a) t a = 0.28, t b = 0.32, t c = 0.4; (b) t a = 0.15, t b = 0.25, 
t c = 0.6. 

associated with the two Dirac points annihilate each other at the critical point. 

3. Entanglement entropy of noninteracting fermions 

The entanglement entropy in the ground state of free fermions can be determined from 
the correlation matrix as a fictitious Hamiltonian defined as [121 Rf3l 05j 

Gij = Tr {pac\cj), (9) 

where q (cj) is the annihilation (creation) operator for site i,j G A. There exist an 
simple relation between correlation matrix of order O(Na) and the reduced density 
matrix of order 0(2 Na ) : pA= det(l — G) exp{^Dfln G(1 — G)]ijc\cj} [®j, where Na is 
the number of sites in A subsystem. Calculating the eigenvalue X n of matrix G requires 
only the diagonalization of Na x Na matrix, and then the entanglement entropy is given 
by 

S vN = —Tr[GTog 2 G - (1 - G) log 2 (l - G)]. (10) 

The eigenvalue X n (0 < X n < 1) of G is termed entanglement spectrum, which is an 
alternative description of Li and Haldane’s version m- Recently the entanglement 
spectrum has been proposed as a ground-state property that captures characteristic edge 
excitations [38l. ; 49l [50] . In a large subsystem, most of the eigenvalues lie exponentially 
close to 0 and to 1 [51 j. A special role is played by eigenvalues X n = 1/2 corresponding 
to a zero-energy mode [52], and it will contribute a maximal value 1 to the entanglement 
entropy by Eq. m- 

In the following, we consider a torus of size L x x L y , and a subregion A can be 
designated as fi= [0, ^iL x ]x[0, z 2 L y \ with z t < 1. Here we measure the entanglement 
between a cornerless cylindrical region A and its complement by two virtual cuts along 
x or y direction; see illustration in Fig[l](a). Zigzag edges emerge following a horizontal 
bipartition, i.e., Z\ — 1, z 2 = 1/2, and k x labels the momentum running along the zigzag 
ribbon. A vertical truncation produces armchair edges, i.e., z\ — 1, z 2 — 1/2, and k y 
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is well defined instead. One can simultaneously classify the entanglement spectrum 
with respect to the momentum parallel to the cut, and an intimate relation between 
entanglement spectrum and eigenspectrum of subblock A can be discovered. 

In Figj3] (a), a flat band of zero energy exists for momenta connecting with two 
Dirac points along x direction near the zigzag edges. Such zero-energy edge states are 
reflected by the maximally entangled modes in entanglement spectrum, as shown in 
Figj3] (b). The fermion doubling in a chiral symmetric system is a universal property 
of the bulk, and the bulk-edge correspondence guarantees the edge states appear as 
pairs. In other words, such doubly degenerate entanglement modes are protected by the 
chiral symmetry. In the semimetallic phase with zigzag edges, the contributions to the 
entanglement entropy can be broadly separated into those from the edges and those from 
bulk, while only bulk states contribute to the entanglement entropy in the insulating 
phase (see Figj4j). However, the situation is quite different for armchair edges. The 
interedge interaction of armchair edges in the semimetallic phase makes the localized 
states fade away [53], and thus edge states play a trivial role in the entanglement entropy 
(see Fig. [5]). Nevertheless, as is shown in Figj6j the anisotropy will induce number of 
edge states, which contribute significantly to the entanglement entropy. It is explicit that 
the zero-energy states have one-to-one correspondences with the maximally entangled 
modes in the entanglement spectrum. In a sense, the entanglement spectrum is a very 
promising tool to characterize chiral edge states. 






k 


X 


Figure 3. (a) The energy spectrum for f a =0.28, fb=0.32, f c =0.40 on L x = 40, L y = 80 
lattice with zigzag edges, (b) The corresponding entanglement spectrum at g = 0. The 
inset amplifies the region in the dashed circle. 


We show the half-space entanglement entropy as a function of t c along the 
normalized path t a + t b + t c = 1 for t a = 0.15, fJL — 0 in FigJTJ The entanglement entropy 
changes discontinuously when t c moves across points t l c = \t a — t b \ and t 2 c = t a + t b . 
In fact, the abrupt change associated with the topological phase transition can be 
grasped more clearly by its first derivative for both types of bipartition and any /j 
within the bulk bands. The system is an insulator when t c < \t a — t b \ or t c >t a + t b , and 
enters the semimetallic phase when \t a — t b \< t c <t a + t b . In the insulating phase, the 
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Figure 4. (a) The energy spectrum for i a =0.15, 0.25, f c =0.6 on L x = 40, L y = 80 

lattice with zigzag edges, (b) The corresponding entanglement spectrum at g — 0. 




Figure 5. (Color online) The energy spectrum for f a =0.28, fb=0.32, f c =0.4 on 
L x = 40, L y = 80 lattice with armchair edges, (b) The corresponding entanglement 
spectrum at g = 0. The inset amplifies the region in the dashed circle. 




Figure 6. The energy spectrum for f a =0.15, tb= 0.25, f c =0.6 on L x = 40, L y = 80 
lattice with armchair edges, (b) The corresponding entanglement spectrum at g = 0. 
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entanglement entropy is linear with the block length for both edges, and the deviation 
seems negligible. However, the fini te-size effect is prominent in the gapless phase. There 
exists a remarkable correction to the linear dependence on the block length. We note 
that the correction is likely a logarithmic form for armchair edges, whereas an oscillating 
correction exists for zigzag edges. We attribute this peculiar correction to finite-size 
effect of edge modes. The trends in the proliferation of edge modes for both edges are 
converse as we approach a quantum critical point, as displayed in the inset of Figd 
We find that the existence of edges state is that t a , h and t c should satisfy the triangle 
inequality for zigzag edge [53] . 



Figure 7. (Color online) The entanglement entropy through path t a + tb + t c = 1 
for t a = 0.15,/u = 0 with zigzag edge (left axis) and armchair edge (right axis). Inset 
shows the convolution of the number of edge states on L x = 30, L y = 60 lattice and 
g = 0 . 

Next we investigate the scaling behavior of bipartite entanglement entropy in the 
presence of different edges. The scaling of the entanglement entropy has attracted much 
attention by a flurry of recent works [55] 561 [57 :, 1.58]. A general form of the entanglement 
entropy in two-dimensional system is argued to be of the form 

S v n = ci L log 2 L + c 2 L + c 3 log 2 L + c 4 , (11) 

where c* (i— 1,- - -, 4) are size independent coefficients. Eq. (JTTj) includes a universal 
logarithmically divergent correction and a nonuniversal area-law term. Which term is 
the leading term requires further elucidation of the Fermi surface. We will give both 
analytical and numerical results for the logarithmical scaling of entanglement entropy 
in one concrete case. 

Apparently in the insulating phase, the entanglement entropy grows asymptotically 
as L [59], as supported in Fig]8](a). For /i = 0 in the semimetallic phase, the Fermi 
surface reduces to zero dimension, and then the entanglement converges to the area law 
[55 , 601 [6TJ . as is shown in Fig]9](a). We see that the slope of scaling changes for varying 
parameters in the semimetallic phase, so the area-law term is nonuniversal. Apart from 
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the leading power-law term, the subleading term is of general interest. To this end, we 
plot the quantity 

= LS vN (L + 2 ) - (L + 2 )S vN (L), ( 12 ) 

which will cancel the area-law term and hence extract the subleading contribution. It is 
remarkable that behave diversely in different phases. If the entanglement entropy 
follows a behavior as 5 v n(T) = c^L + C 4 + 0(1/L), 5’*^’ will approach a constant term 

like 


S3? ~ -C4, (13) 

as is shown in the inset of FigJHl(a). In contrast, SVn displays a divergent oscillation in 
the inset of Figj9] (a). A careful analysis identifies a logarithmic additive term in these 
critical phases by following S V ^(L) = C 2 L + C 3 log 2 L + 0(1), and then the subarea law 
of the entanglement entropy yields[61, [62] 

SvN ~ —c 3 l°g 2 L. (14) 


The subdominant contribution stems from one chiral mode in the zero-dimensional Fermi 
surface, which infers C 3 = 1/3 by the FisherCHartwig conjecture [631164] . The oscillation 
can be interpreted as the following: I 11 case of symmetric hoppings, the existence of zero- 
energy localized states is restricted for 1/3 of the total momentum parallel to the zigzag 
edge, and hence a multiple of 3 in block length favors more edge states. Hence, the 
oscillation has a period of 3, as is observed in the inset of FigJT] 

A significant difference of the scaling in the entanglement entropy is found in the 
presence of finite density of state at the Fermi surface. I 11 such situation, every low- 
energy radial mode is approximately treated as a chiral excitation that contribute a I 11 L 
to the total entanglement entropy [65] • Therefore, a multiplicative logarithmic correction 
to the power-law behavior arises for finite Fermi surface (i.e., line nodes), which are 
reported in Fig] 8 ](b) and Fig]9](b) [ 66 , [67] . Noticeably, we see that the entanglement 
entropy with armchair edge is larger than that with zigzag edge, and the reason will be 
accounted soon. 

Widom conjecture [ 68 ] predicted that the prefactor of the logarithmic term is 
determined by chemical potential p via an accurate formula [69] 


ci(h) = 


24t r 


dS T 


dS k |n x • n k | 


ion 


'MV) 


(15) 


where H is the real-space region of A, rescaled by L such that Vol(f2)=l. Vectors n x 
and n k denote the normal vectors on the real-space surface dfl and the Fermi surface 
<9r(/i). From Eq. (fT5lh we see the topology of the subsystem and the Fermi surface 
has intimately related to the entanglement entropy. Thus the existence of edge states 
crucially affects the physical properties of graphene system. Besides, it was argued 
that Lifshitz quantum phase transitions are accompanied with a nonanalyticity in the 
prefactor of the leading-order term (SI- In two dimensions the calculation of c 1 (/ 1 ) can be 
simplified. The normal unit vectors in real space are ±y for zigzag-type bipartition and 
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Figure 8. (Color online) The scaling of entanglement entropy at t a =0.15, tb=0.25, 
f c =0.60 when (a) g = 0 and (b) g = —0.5. Lines are linear fits. for f o =0.15, 

tb=0.25, i c =0.60 with zigzag edge is shown in the inset. 




Figure 9. (Color online) (a) The scaling of entanglement entropy with zigzag (square) 
and armchair edges (circle) at t a =tb=t c =1.0 (solid) and f a =0.28, tf,=0.32, t c =0.40 
(open) with g = 0. for t a =t b =t c =1 with zigzag edge is shown in the inset, (b) The 
scaling of entanglement entropy with both zigzag and armchair edges at f o =fb=t c =1.0 
with g = —1.0. Lines are linear fits. 


±x for armchair-type bipartition. Since the Fermi surface is a curve in two dimensions, 
we can parameterize the curve by k x = k x {9) and k y = k y {6) like [Mj 


(t a + t b ) cos k x + t c cos k y = /i cos 9 , 
(t a — t b ) sin k x — t c sin k y = p sin 9. 


Then we have 


c i,zigzag (f) 

c i, armchair (f) — 


1 [ , 

dk x 

- / d9 


J-27T Jon 

d9 

1 f 

dk n 

/ d9 


J-2vr J d n 

d9 


(16) 

(17) 


Eq. (IT6l) and Eq. ([T7|) indicates Ci(/i) are proportional to the variation of Fermi surface 
along x and y , respectively. For \p\ is large, the Fermi surface is an oval, while 
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the disconnected Fermi components appear with the decrease of \p\. Without loss 
of generality, we give a simple demonstration of the solution for symmetric case 
t a = h = t c = 1. The dispersion has particle-hole symmetry and 90-degree rotational 
invariance, so we only consider positive p in the first quadrant of the Brillouin zone, 
and then the Fermi surface can be parameterized by 

k x (9) = cos -1 R(p, 9) 
k y {6 ) = sin _1 (— psm.9) 

for 0 < p < 1 and 6 G (—vr,0). Here R(p,9) = (— \/l — p 2 sin 2 6 + pcos9)/2. 
The longitudinal variation A k y is 2 arcsin(yu) and the horizontal variation A k x is 
arccos( "^ 1 ) — arccos (^-). When p > \/5, it is replaced by 

k x {0) = cos -1 R(p, 9) 
k y {9 ) = sin -1 (—p sin 9) 


for 9 G [— arccos(^^), 0]. The longitudinal variation A k y is arccos 
horizontal variation A k x is arccos (^-). Otherwise when 1 < p < \/5, 

k x (9) = cos -1 R(p, 9) 
k y (9) = sin _1 (— psin9) 


F 2 - 5 


and the 


for 9 G [— sin 1 -, 0] and 


k x {9) = 7r — cos 1 R(p, 9) 
k y (9 ) = 7T + sin -1 (p sin 9) 

for 9 G [arccos— 7T, arcsin(^) — 7r]. The longitudinal variation A k y is arccos iL ^ 1 
and the horizontal variation A k x is arccos( i£ 2 ^). The resulting prefactors of leading 
logarithmic terms are: 


'i,zigzag 


(aO q ^ k x , 

37r 


c i, armchair (aO 


= —A k v . 
3tt v 


(18) 


A generalization to anisotropic case and arbitrary edges is straightforward. The 
coefficient of the leading term in system of open boundary condition is half of the one of 
periodic boundary condition, due to the fact that the boundary between the bipartite 
spaces is half of that of stripes m 

From Fig. [TO], we can recognize the entanglement entropy presents a shoulder-like 
behavior with respect to the chemical potential p. The finite-size effect becomes less 
dominated for large size. Ci(/i) extracted from fits to our numerical data agrees rather 
well with the exact form Eq. (fT8|) in the inset of Fig[10] C\{p) goes to zero with p 
approaches Dirac point, where the logarithmic correction is vanishing. When p deviates 
slightly from the critical point, two small Fermi pockets are formed with Fermi vector 
kp = \p\/vf, and Eq. dlSli implies that Ci(/i) scales linearly for the low-energy modes 
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Figure 10. (Color online) The entanglement entropy between A and B block for 
different chemical potential at t a =tb=t c = 1.0 with (a) zigzag edges and (b) armchair 
edges. The legends mark the half size. Insets show the prefactor c\{g), 02 (g) in Eg. (fill) 
as a function of the chemical potential g for the ground state of the two-dimensional 
fermionic tight-binding model. c\(g) extracted from fits to our numerical data (red 
filled square) are compared with the analytical results (solid line) in the upper left 
inset, (c) Analytical results c\(g) for both edges, and the ratio in the inset. 


when system size L l*, where l* is a characteristic length l* ~ kf . The ratios 
of c\(n) to C 2 (/i) with both edges are above a factor \/3, because the Brillouin zone 
of brick-type lattice along y direction is stretched \/3 times longer than that along x 
direction comparing with pristine honeycomb lattice. The ratios reach their maxima 
passing through a Van Hove singularity, i.e., fi = 1, where two Fermi surfaces merge 
to form a single Fermi surface; see illustration in Fig JTOf c). Thus at /i = 1 there is a 
Lifshitz transition. Besides a logarithmic divergence of the entanglement entropy, the 
linear-term coefficient is intricately related to ultraviolet cutoff or correlation length, 
and it is generally nonuniversal. Shown in the lower insets of FigJTO] (a-b) are C 2 (/i), 
which decrease with the increase of |/r|. Extraordinarily, there is a dip at p = 1 by a 
meticulous observation. We attribute such decrease to the enhanced contribution from 
a negative subarea term at /i = 1 [61] . 

4. Spin-orbit interaction 

Next, we extend our scope to Kane-Mele model, which is a tight-binding model with 
second neighbor spin-orbital interaction defined on a honeycomb lattice, 

H — — t "y ^ y ^ { a \cPr+Ti,a + b r +Ti,a a r,a) 

rEA^a i=a,b,c 

+ it 2 E E a r,a[ Cr ' * '^j)]a,/3 a r+n+Tj,l3 

rGA^4,a,/3 i,j=a,b,c 

- it 2 b rA a ' x 

rEA B-,ot,f3 i,j=a,b,c 
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t 1 ^ 1 a r,a a r,a h ^ ^ ^r,oPr,a- (19) 

rEAs,a 

Here f\ = —r* is the antiparallcl vector of T; in Eq.([4]). The Hamiltonian is time 
reversal invariant and has a wide range of applicability, such as quantum spin-Hall 
(QSH) insulator [72] [73], graphene/graphane interface [74] and silicene [75] 7611771 ^78]. 
The energy band structure of the brick lattice is given by 

<* = ±7*2 + e>2, (20) 

where 4> fc = —t(e lky + e lkx + e ~ lkx ) and 0 fc = 4f 2 cos k y sin k x . Without the second- 
neighbor hopping, the model’s bulk and surface spectra display the non-generic feature 
that all energy levels come in pairs at each momentum due to parity symmetry [52]. 
The presence of the spin-orbit interaction, i.e., f 2 7^ 0, leads to a gap in the bulk 
spectrum, as shown in FigJTTl More interestingly, in the case of zigzag ribbon, the flat 
dispersion is replaced by the helical modes intrinsic to the QSH insulator as the spin- 
orbital is introduced im The dominant contributions originates from the helical edge 
state is displayed in FigJTTTd). The area-law monotonic scaling at p = 0 is confirmed 
in Fig JTTT a). and the leading logarithmic divergence is also observed in Fig JTTT b). 





Figure 11. (Color online) (Color online) The energy spectrum of Hamiltonian (TTH1) 
for (a) t 2 = 0 and (b) f 2 = 0.1 on L x = 40, L y = 40 lattice with zigzag edges. The 
corresponding entanglement spectrum at y = 0 for (c) f 2 = 0 and (d) f 2 = 0.1. 
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Figure 12. (Color online) The scaling of entanglement entropy with zigzag edges at 
t2=0 and 0.1 when (a) g = 0 and (b) g = —1. Lines are linear fits. 


5. Conclusion 

In this paper, we study the von Neumann entropy in free fermionic systems on a two- 
dimensional honeycomb lattice. Distinct bipartitions on the honeycomb lattice will 
produce different edges, and the edge states will have an effective impact upon the area- 
law terms in the entanglement entropy. We exhibit that the zero-energy edge state has 
an one-to-one correspondence with the maximal entangled mode in the entanglement 
spectrum, and hence we speculate the main contributions to the entanglement entropy 
come from the bulk states and edge states. The edge states play an essential role 
in semimetallic state with zigzag edges and in insulating phase with armchair edges. 
Furthermore, we find the entanglement entropy obeys area law when the Fermi surface is 
within the gap or a nodal point. The corrections to the pow-law term behave differently 
in critical phase from in noncritical phase. A logarithmic violation to the area law 
emerges in the presence of line-nodes Fermi surface. We also show that the logarithmic 
scaling highly depends on the topology of the Fermi surface. The prefactor of logarithmic 
term is determined by the variation of Fermi surface along the direction parallel the edge. 
The scaling laws are independent of the statistics of the microscopical constituents. The 
introduction of spin-orbit coupling will cause a gap in the bulk spectrum and a decrease 
of entanglement entropy. 

Since there is an equivalence between the entanglement entropy and particle number 
fluctuations in free-fcrmion systems pl Ea EH], it is very promising to measure the 
entanglement entropy in the graphene-like systems. Especially experimentally the edges 
can be better fabricated in the artificial honeycomb than those of natural electronic 
graphene, which tend to be very irregular and contaminated with adsorbates. The 
flat edge modes can be observed in the zero-bias conductance measurements [80]. Our 
results can be extended to other free fermionic systems, i.e., band topological insulators 
and superconductors, as well as weakly interacting systems [81 ]. 
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